function IQ=pwdoppler%(prf,v,fs,fc,B);
% Pulsed wave signal from moving reflector
z0=0.01;%start position
z1=0.013;%end position
v=0.5;%m/sec
prf=20e3;%pulse repetition frequency

complexDemod=1;%=1 to perform complex demodulation

prt=1/prf;
dz=v*prt;

zs=0.0125;
z=z0;
figure(1);clf;
pulseEchoPlot(z,zs,complexDemod);
pause;
n=1;
while z<z1,
    figure(1);
    [IQ(n,:),iq(n)]=pulseEchoPlot(z,zs,complexDemod);
     t(n)=n*prt;
    %figure(2);plot(iq(n),'*');axis([-1,1,-1,1]*1.5);grid;pause(0.1);
    n=n+1;
    z=z+dz;
end;
figure(2);clf;
if complexDemod,
    subplot(2,1,1);plot(iq);axis('image');axis([-1,1,-1,1]*1.5);grid;
    xlabel('Real part');ylabel('Imaginary part');
    subplot(2,1,2);plot(t,real(iq),t,imag(iq));
else plot(t,real(iq));
end;

function [iq,iq1]=pulseEchoPlot(z1,zs,complexDemod);
%Pulse echo 1D imaging - simulation
% z1 is distance to reflector [m]. default=0.01 m
%23.01.01  Hans Torp
% Modified for complex demodulation 17.10.01

if nargin<1,z1=10e-3;end;

c=1540;%speed of sound
fs=40e6;%sampling freq.
tmax=20e-6;
t=0:1/fs:tmax;%time axis
z=c/2*t;dz=z(2)-z(1);%depth axis

%define  transmit pulse
f0=5e6; B=0.6e6;
Tp=2e-6;%pulse length
tp=0:1/fs:Tp;
p=sin(2*pi*f0*tp);%rectangular pulse
p=XmitPulse(fs,f0,B);
obj=0*z;obj(1)=2;%object function;
obj(1+round(z1/dz))=1;%reflector in position z1
s=filter(p,1,obj);%received signal

figno=1;
figure(figno);set(figno,'DoubleBuffer','on');%prevent flickering

sp=hilbert(s);%pre envelope
sh=imag(hilbert(s));%hilbert transformed signal
smix=exp(-i*2*pi*f0*t);%down-mixing signal
iq=sp.*smix;%in phase and quadrature components
[dum,is]=min(abs(z-zs));
iq1=iq(is);
ts=2/c*zs;

myaxis=[t(1),t(end),-2.2,2.2];
subplot(4,1,1);plot(z,obj);xlabel('depth z [m]');axis('tight');ylabel('Object');
subplot(4,1,2);plot(t,s);xlabel('time [s]');axis(myaxis);ylabel('RF-signal');
hold on; plot([ts,ts],[-1,1]);hold off;
if complexDemod,
    subplot(4,1,3);plot(t,sh);xlabel('time [s]');axis(myaxis);ylabel('Hilbert transf.');
    hold on; plot([ts,ts],[-1,1]);hold off;
    subplot(4,1,4);plot(t,real(iq),t,imag(iq));xlabel('time [s]');axis(myaxis);ylabel('IQ signal');
    hold on; plot([ts,ts],[-1,1]);hold off;
end;

pause(0.01);

function pxd=XmitPulse(fs,fc,B,Tp);
%define transmitted pulse
%fs: sampling frequency
%fc: pulse center frequency
%B  pulse bandwidth
% Tp  pulse length, default Tp=3/B;
if nargin<4, Tp=3/B;end;
[bxd,axd]=butter(2,2*[fc-B/2,fc+B/2]/fs);%butterworth bandpass filter
dt=1/fs;
txd=0:dt:3e-6;txd=txd';
impuls=zeros(size(txd));impuls(1)=1; %
pxd=filter(bxd,axd,impuls);
pxd=pxd/max(pxd);


